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ABSTRACT 

We present a statistical exploration of the parameter space of the De Lucia and 
Blaizot version of the Munich semi-analytic (SA) model built upon the Millennium 
dark matter simulation. This is achieved by applying a Monte Carlo Markov Chain 
method to constrain the 6 free parameters that define the stellar and black-hole mass 
functions at redshift zero. The model is tested against three different observational 
data sets, including the galaxy /\-band luminosity function, B — V colours, and the 
black hole-bulge mass relation, separately and combined, to obtain mean values, con- 
fidence limits and likelihood contours for the best fit model. Using each observational 
data set independently, we discuss how the SA model parameters affect each galaxy 
property and find that there are strong correlations between them. We analyse to what 
extent these are simply reflections of the observational constraints, or whether they 
can lead to improved understandings of the physics of galaxy formation. 

When all the observations are combined, we find reasonable agreement between 
the majority of the previously published parameter values and our confidence limits. 
However, the need to suppress dwarf galaxy formation requires the strength of the 
supernova feedback to be significantly higher in our best-fit solution than in previous 
work. To balance this, we require the feedback to become ineffective in halos of lower 
mass than before, so as to permit the formation of sufficient high-luminosity galaxies: 
unfortunately, this leads to an excess of galaxies around L* . Although the best-fit is 
formally consistent with the data, there is no region of parameter space that reproduces 
the shape of galaxy luminosity function across the whole magnitude-range. 

For our best fit we present the model predictions for the 6j-band luminosity and 
stellar mass functions. We find a systematic disagreement between the observed mass 
function and the predictions from the A"-band constraint, which we explain in light 
of recent works that suggest uncertainties of up to 0.3 dex in the mass determination 
from stellar population synthesis models. 

We discuss modifications to the semi-analytic model that might simultaneously 
improve the fit to the observed mass function and reduce the reliance on excessive 
supernova feedback in small halos. 

Key words: methods: numerical - methods: statistical - galaxies: formation - galax- 
ies: evolution 



1 INTRODUCTION 

The combination of the well established ACDM paradigm 
and the geometrical growth of computer power in recent 
years, has allowed direct N-body simulations to give us a 
comprehensive picture of the formation and evolution of 
dark matter structure, from the primordial gravitational in- 
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stabilities to the form ation of massive sup erclusters today 
(The Millennium Run. ISpringel et alj|2005h . 

However this progress is not reflected in our theoret- 
ical understanding of the behaviour of the baryons, with 
much of the physics governing galaxy formation and evolu- 
tion still poorly understood. A combined simulation of gas 
and dark matter resolving small-scale galaxy evolution pro- 
cesses (such as gas cooling, star formation and feedback) 
over a cosmologically interesting volume is still years away. 
Therefore the only plausible method to construct a model 
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galaxy population for comparison with observed large-scale 
surveys is us ing a semi-ana l ytic f ormalism. Originally in- 
troduced by IWhite fc Reesl (1 19781 ). this formalism treats 
the dark matter structu re using either a Press-Schechter 
ft Press fc Schechtedll974h . Monte-Carlo or N-body approach 
and on top of that follows galaxy evolution using param- 
eterised equations governing the laws of subgrid physics. 
The b asic methodology was set bylColel (|l99ll ). lLacev fc Silkl 
i|l99ll ) and lWhite fc Frenkl (|l99lf ). including the dependence 
of gas cooling and star formation on the dark matter halo 
density profile, feedback and chemical enrichment to account 
for the effect of supernova (SN) explosions on the properties 
of the hot gas, and stellar population synthesis models to 
convert star formation histories into observed stellar prop- 
erties. 

Further developments included newly-derived stellar 
population mo dels and improved star formation and SN 
feedback laws (Kauffmann et all 1 19931 ; lLacev et alj 1 19931 ; 
ICole et al.l 1 1994) . W th the level of complexity achieved, 
they were able to predict a large range of galaxy properties 
such as star-formation rates, luminosity functions and rela- 
tions between circular velocity, luminosity, metallicity and 
mass-to-light ratios. 

From this original recipe, two models started evolv- 
ing separately, one mainly based in Munich and another 
in Durham. By the end of the decade, m ost of the modern 
day prescriptions were already introduced (jKauffmann et all 
I1999) ; ICole et all l2000h , including gas cooling, star forma- 
tion, chemical enrichment and dust extinction, calculations 
of disks and bulge properties, stellar population synthesis 
models, merger follow up with dynamical friction and an 
early ve rsion of the SN feedback treatm ent. 

In iKauffmann fc Haehneltl <|200d ) a model for the 
growth of black holes (BHs) d ue to instabili t ies a rising 
from mergers was pro posed, and iBenson et al.l (|2003l ) and 
iDe Lucia et all Swoj) studied new treatments of the SN 
feedback, including the current model where the SN can 
not only reheat the cold gas into the hot phase but also 
eject gas from the halo (to be reincorporated at later 
times). In parallel, a number of other groups begun to 
develop independ ent models, to study different aspects of 
galaxy formation dSomerville fc Primac k 1999; Menci et al.l 



20021; lHatton et al.ll2003l ; iDaigne et al.ll2004 ; iMonacd 12004 
Kane et al.ll2005l ). 

Most recently, the Munich and Durham semi-analytics 
have been combined with the Millennium dark matter sim- 
ulation and an additional recipe, the BH radio mode, in- 
troduced to reproduce the quenching of gas cooling star 
for mation in the gas surrounding central c l uster galax- 
ies ilSpringel et aljBool ICroton et all 120061; iBower et all 



20061; see al s o ICranato et all 12004 ICattaneo et all l2006t 
Menci et all 120061 ; iMonaco et all 120071 ; ISomerville et all 
200Sh . 

Finally, present day studies include new dust models 
|De Lucia fc Bl aizot 2007J), the st udy of alternative feed back 
processes such as galactic winds (|Bertone et alj|200^ ). im- 
proved recipes for the stripping of gas during galaxy mergers 
l|Font et al.ll2008h and investigation of the ability of the en- 
ergy released by AGN feedb ack to reproduce th e properties 
of the intra-cluster medium |Bower et al.ll2008l ). 

With all these recipes in place, the models success- 
fully reproduce a vast range of observable properties, from 



galaxy luminosities and colours, including environment de- 
pendences, to scaling relations such as Tully-Fisher dia- 
grams. However, until now the level of agreement with ob- 
servations and the relative weight of different observations 
in the final choice of the parameters has never been studied 
in a statistically-consistent way. 

Moreover, the large number of observational properties 
that the models aim to predict requires a large number of 
parameters (some of which are strongly correlated) , produc- 
ing considerable difficulties in determining how to improve 
the agreement with new observations without destroying the 
match with existing data sets. In addition, whenever reason- 
able agreement proves to be impossible, it is hard to know 
whether there is a failure in determining the right parameter 
configuration, whether there is a fundamental problem with 
the underlying model, or whether the introduction of new 
physics is called for. 

These difficulties can be overcome by combining multi- 
ple observations with proper sampling of high-dimensional 
parameter spaces. This has proved to be a fruitful approach 
in theoretical cosmology where techniques such as Monte 
Carlo Markov Chain (MC MC) paramet er estimation have 
been extensively used (see iTrottal (|2008l ) for a comprehen- 
sive review). The aim of this paper is to introduce MCMC 
techniques into semi-analytic models of galaxy formation. 
While thi s work was being dev e loped , a first result was pro- 
duced by iKampakoglou et all l|2008h . These authors have 
introduced simi lar tools to their ow n semi-analytic recipe, 
an extension of IDaigne et all (12004 ). which uses a statisti- 
cal method to generate halos. We differ from them in tha t 
we use the semi-analytic model IDe Lucia fc Blaizotl l|2007h . 
in our case built upon a direct dark matter simulation of a 
cosmological size (the Millennium Run). 

This will allows us to understand how galaxy proper- 
ties are affected by individual parameters, obtain confidence 
limits for the parameters, and verify the agreement between 
the model and different observations in a statistically robust 
way. 

This paper is organised as follows. In section[2]we briefly 
describe the semi-analytic model used in our study. In sec- 
tion [3] we present the MCMC technique used to constrain 
the model parameters and explain how it is implemented 
into the semi-analytic recipe. Section [4] describes the obser- 
vational data used in this work, clarifying which parameters 
are constrained by each observational galaxy property. In 
section [5] we present our results including correlations be- 
tween the parameters analysed and predictions for our best 
fit model. Finally in section[6]we summarise our conclusions. 



2 THE MODEL 

In this section we brie fly describe the sem i -analy tic model 
we use for this work |De~ Lucia fc Blaizotl 120071 hereafter 
DLB07), and the unde rlying dark matter simulation, the 
Millennium Simulation (|Springel et al.ll2005h . 

The Millennium Simulation traces the evolution of dark 
matter haloes in a cubic box of 500/t _1 Mpc on a side. It 
assumes a ACDM cosmology with parameters S7 m = 0.25, 
fib = 0.045, h = 0.73, Ov = 0.75, n = 1, and a 8 = 0.9, 
where the Hubble parameter is Ho = 100 hT 1 kms -1 Mpc -1 . 
The simulation follows 2160 3 dark matter particles of mass 
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Table 1. Best-fit parameters for the SA model from DLB07. The 
first 6 parameters are frozen in our analysis at the values shown 
here. A detailed description of the parameters is given in the text. 
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8.6 x 10 h~ Mq. Since dark matter haloes are required 
to contain at least 20 particles, the minimum halo mass is 

1.7 x 10 10 h^ 1 Mq, with a corresponding baryonic mass of 
about 3.1 x 10 9 /i _1 M©. 

The coupling of the semi-analytic model onto the high- 
resolution N-body simulation f ollows the technique imple- 
mented bv lSpringel et all |200lf) . The treatment of physical 
processes driving g alaxy evolution builds o n the methodol- 
ogy in trodu ced by Kauffmann et al.| (| 1999T ) , ISpringel et al.l 
(200 1| ) and |Pe Lucia et al l (|2004ll. and is a s li ghtly mod- 
ified version of th at used in ISpringel et all i|2005t l and 
ICroton et all (|2006l ). 

The model divides the gas content of galaxies into sev- 
eral distinct phases. When the galaxy first forms, the gas 
enters the hot halo at the virial temperature. It can then 
cool down to join the cold disk. Stars form from the disk 
and feedback of energy via supernovae can cause gas to heat 
back up from the cold disk to the hot halo. Finally, it is 
possible to heat gas still further, ejecting it from the galaxy 
into an external reservoir from which it gradually leaks back 
into the hot halo. 

From the original 12 parameters in the model, we choose 
to freeze 6 of them at the values chosen by DLB07, as shown 
in the top row of table [11 The cosmic baryon fraction, 
is fixed by the cosmologjij, while the redshifts of the begin- 
ning and end of reionization (zo, z r ) are used to modify the 
baryon fraction in small halos, accounting for the effects of 
phot o-ionizing heating dKravtsov et al ] |2004 I Croton et al.l 
|2006| ). Tmergcr is the threshold mass ratio that defines the 
distinction between Major and Minor mergers. R is the re- 
cycled fraction and Y the yield, both of which depend upon 
the details of the stellar initial mass function. The obser- 
vational data that we use in this work do not allow us to 
strongly constraint any of these values. 

That leaves 6 free parameters in our study. These are 
the star formation efficiency, o?sf, the fraction of cold gas 
accreted by the central BH during mergers, frm, the quies- 
cent hot gas BH accretion rate, kAGN, the SN feedback disk 
reheating efficiency, eaisk, the SN feedback halo ejection effi- 
ciency, ehaio, and the ejected gas reincorporation efficiency, 
7ej. We briefly describe the meaning of each of these p aram - 
eters below, for a full description see ICroton et alj (|2006l ) 
and DLB07. 

The model converts cold gas into stars at a rate given 



1 Note that, as in DLB07, we use the value of 0.17 suggested by 
WMAP rather than the value of 0.18 used to generate the power 
spectrum for the Millennium Simulation 
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^dyn.disk 
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where m co id is the mass of cold gas, m cr it is the mass that 
corresponds to a critical surface den sity above whic h gas 
can collapse and form stars (following lKennicuttlll998h , and 
£dyn,disk is the dynamical time of the disk. Note that the 
fraction of mass locked up in stars is (1 — R)rh- tc dt, the rest 
being instantaneously returned to the disk. 

As massive stars complete their life cycle, SN events 
start injecting energy into the surrounding medium, reheat- 
ing cold disk gas and even ejecting gas from the hot halo. 

For each mass Am* turned into stars, the amount of 
gas reheated from the cold disk to the hot halo is given by 



Am re heated — EdiskAm* 



(2) 



with the cano nical efficienc y of 3.5 being motivated by ob- 
servations bv lMartinl (|l999h . 

This proves insufficient to prevent star-formation in 
dwarf galaxies as the cooling times are so short that the 
gas rapidly cools back down to rejoin the disk. For this rea- 
son, and moti vated also by observations of galactic outflows 
l|Martinlll996T ), the models allow SN to expel gas completely 
from low-mass galaxies. 

The amount of energy released by supernova during the 
formation of Am* stars is 

A£ S n = 0.5 ehaio AjtuVsn (3) 

where Vsn = 630 km s" 1 . Any escesfl energy left over from 
reheating the cold gas is used to eject a mass of gas Am e j e cted 
from the galaxy 
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where V wiT is the circular velocity of the dark matter halo. 

This ejected gas is kept in an external reservoir and 
returned to the hot halo at a rate 

^'ejected 



^/ejected 



-7ej" 



dyn 



(5) 



where rejected is the mass of ejected gas and tdyn is the 
dynamical time of the halo. 

SN feedback is ineffective in large galaxies and so an- 
other form of heating must be included. Without it, cen- 
tral cluster galaxies appear too massive and too blue, an 
aspect of the cooling flow problem, well-known to X-ray as- 
tronomers. The solution is thought to be mechanical heating 
by BHs accreting at well below the Eddington limit, with the 
amount of energy released depending on the mass accretion 
rate of the central supermassive BH, which in turn depends 
on the BH mass. 

To describe this process it is necessary to introduce 
two different modes of AGN activity: the quasar and radio 
modes. The former is thought to be inefficient at heating the 
gas but is primarily responsible for BH growth. It was origi- 
nally introduced into the models simply to predict the mass 
of central BHs and the corresponding SA model parameter, 



2 Note that the reheated fraction is not reduced if this excess is 
negative. 
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/bh, regulates the BH growth by accretion associated with 
galaxy mergers 

a /BH(m sa t/m C entr a l) m C old , r .s 

AmeH.Q = 1 + (28Qkms - 1/Kir)2 ■ (6) 

The radio mode reflects the BH growth via quiescent ac- 
cretion in a static hot halo. It may represent either Bondi ac- 
cretion directly from the hot phase, or the accretion of small 
quantities of cold gas. It is described by the phenomenolog- 
ical model 

3 



772BH,R = ^AGt 
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where mBH is the BH mass and /hot is the mass fraction of 
hot gas in the halo. 

The radio mode makes a minor contribution to the 
growth in mass of the BH but is assumed to generate me- 
chanical heating at a rate 



Lbh = r]mBH,R c 



(8) 



where c is the speed of light and the efficiency parameter r\ 
is frozen at 0.1 □ This heating is used to reduce the rate at 
which gas cools from the hot halo into the cold disk. 



3 MONTE CARLO MARKOV CHAIN 

3.1 Bayesian Monte Carlo Markov Chain Analysis 

Monte Carlo Markov Chain (MCMC) methods are a class 
of algorithms for sampling a multidimensional space with 
a probability proportional to the likelihood that the model 
describes the observationa l constraints . The following brief 
description follows that in IPress et all i|2007l ,. 

A typical application of this method is when it is possi- 
ble to calculate the probability, P(D\x), of a given data set, 
D, given the values of some model parameters, x. Bayes' 
theorem says that, given a prior P{x), the (posterior) prob- 
ability of the model (which will be sampled by the MCMC) is 
ir(x) oc P(D\x) P(x) with an unknown normalising constant. 
The advantages of MCMC are that the posterior distribution 
and correlations for the parameters in study can be easily 
recovered from the sample list and the un-normalised prob- 
ability, and that the computational power required scales 
only linearly with the number of parameters. 

The MCMC method uses a Markov chain to step from 
one point in the sample space to the next, meaning that each 
point is chosen from a distribution that depends only on 
the preceding point (the ergodic property). The transition 
probability p(x2\xi) for stepping from point xi to point xi 
should satisfy the detailed balance equation, 

7r(a;i)p(a;2|xi) = ir(x2)p{xi\x2) (9) 



3.2 Metropolis-Hastings Algorithm 

There are several algorithms that can produce a chain 
with the required properties, the most common being 
the Metropolis-Hastings algorithm l|Metropolis et af1ll953l ; 

3 Note that it is the combination 77 &agn that determines the 
heating rate so that the value of r\ is unimportant. 



a(xi,x 2 c) = min 1 



lHastingdfilfflih . This method requires a proposal distribu- 
tion q(x2\xi) that can assume various shapes, as long as the 
chain reaches everywhere in the region of interest. However, 
an inappropriate choice can delay significantly the conver- 
gence of the chain. Considering the underlying probability 
distribution of our parameters we choose a log-normal pro- 
posal distribution with a width that assures that the final 
acceptance of the chain is between 10% and 40%. 

The chain is then started in a randomly selected point in 
parameter space xi. A new candidate point X2c is selected by 
drawing from the proposal distribution, and the acceptance 
probability a(xi,X2c) calculated using the formula, 

K{X2 c ) q ( Xl \x2 c )\ 

7r(a;i)5(a;2 C |xi) / 

The candidate point is accepted with probability a(xi,X2c) 
and X2 is set equal to X2 C , or rejected and the point left 
unchanged (X2 = Xi). 

The ratio q(xi\x2c)/q(x2c\xi) in equation (|10|l repre- 
sents the prior which we assume to be log-normal. 

3.3 MCMC applied to the Semi- Analytic Model 

Implementing the MCMC sampling approach on the semi- 
analytic model parameter space raises considerable issues 
related not only to the copious amount of I/O (the original 
recipe reads in the full Millennium dark matter trees), but 
also to the volume of calculations required to follow the evo- 
lution of over 20 million galaxies in a cosmological volume, 
through more than fifty redshift slices. At each MCMC step 
the semi-analytic model needs to be run with the proposed 
set of parameters, to compute the acceptance probability 
by comparing the outputted galaxy properties with the ob- 
servational constraints. The size of the calculation and the 
number of steps required for convergence makes it unfeasible 
to perform our analysis using the full Millennium volume. 

The structure of the Millennium Simulation provides a 
way to circumvent this difficulty. The output is divided into 
512 files which have self-contained trees, with the galaxies 
on each treated independently. We choose to perform our 
analysis in a single file with a mean density and luminosity 
function analogous to the full Millennium box. This assures 
that the parameter study done on it is representative of the 
full data set. Only the largest galaxies with stellar masses 
greater than about 10 h Mq are not properly sampled 
this way. 

For our best-fit parameters, we rerun the SA model on 
the full simulation: these results are presented in section [S] 
below. 



4 INDIVIDUAL OBSERVATIONAL 
CONSTRAINTS 

4.1 Overview 

The traditional semi-analytic approach, is to adjust param- 
eters only considering observations at redshift zero. Follow- 
ing this philosophy we select 3 independent and local ob- 
servational data sets: the K-h&nd luminosity function, the 
colour-stellar mass relation and the BH-bulge mass relation, 
to fully constraint the 6 parameters defining galaxy masses 
and formation rates of stars and AGN. 
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In this section we present the observations used in our 
analysis and we show how each individual property con- 
straints the different parameters by running the MCMC 
sampling with one observational data set at a time. The 
output is analyze d using getdist, part o f the COSMOMC 
software package (jLewis fc Bridle! l2003 ). adapted to pro- 
duce Id and 2d maximum likelihood (profile) and MCMC 
marginalised (posterior) distributions. For the independent 
observational properties we use different statistical tests 
to assess the likelihood of the model, reflecting the obser- 
vational uncertainty and the nature of the relation under 
study. 

Running our sampling technique with separate obser- 
vational data sets one at a time allows us to gain insight 
on which degeneracies between the different parameters are 
broken by each additional observation. We start by studying 
the influence of varying the parameters on the final A"-band 
luminosity function. 
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4.2 The if-band Luminosity Function 

Despite being one of the most fundamental properties of a 
galaxy, stellar mass is not easily derived from observations. 
To get estimates for this quantity from the observed lumi- 
nosities it is necessary to assume mass-to-light (M-L) ratios 
based on stellar population synthesis models that include 
still poorly-understood dust corrections, initial mass func- 
tions (IMFs) and metallicity evolution. On the other hand, 
semi-analytic models directly predict mass, but to produce 
observable luminosities the same crudely-established process 
must be taken in the reverse direction. 

This difficulty leads us to use the isT-band luminosity 
function. The _K"-band is known as a good mass indicator 
as it is relatively unaffected by dust and represents a fair 
sample of the stellar population. We combine three observa- 
tional studies i|Cole et al.ll200ll ; iBell et al.ll2003l ; Ijones et all 
l200fj ), respectively from 2DFGRS, 2MASS and 6DFGRS, 
from which we build a final luminosity function. The final 
data points are given by the average of the maximum and 
minimum number-density estimates in each magnitude bin, 
with errors Oi equal to half the difference between them. 

The comparison between the A"-band luminosity func- 
tion from the original DLB07 model (using the published 
parameter values) with the observations is shown in FigfT] 
The original model already shows good agreement with the 
combined data except for the faint end, over predicting the 
number of dwarfs galaxies with magnitudes fainter than 
K » -22. 

To compute the likelihood of the model for the A"-band 
luminosity function we use the chi-square probability func- 
tion where 



X 



(Ni 



Tli) 



(11) 



is summed over the observational bin range plotted in Fig[T] 
and Ni and m represent the number of observational and 
simulated galaxies in each bin, respectively. 

In Fig[2] we plot the la and 2a preferred values from 
the MCMC (solid lines) and the maximum likelihood value 
sampled in each bin (colour contours), for the subset of the 
original parameters (with values plotted in log space) con- 



Figure 1. The galaxy A'-band luminosity function at z=0 for the 
DLB07 model (dashed red line) and our best fit model (solid red 
line) . The model predictions are compared with obse rvations from 
jCole et al.ll200ll : IBell et al.ll2003l : Ijones et af1l2006l) combined to 
produce a new luminosity function reflecting the scatter between 
them. 



strained only by the observational K-band luminosity func- 
tion. 

In interpreting this and future plots, one should bear 
in mind that the contours follow the MCMC sampling in 
parameter space, which should trace out the relative likeli- 
hoods of different regions (the posterior distribution). The 
colours represent the maximum likelihood projected along 
all the hidden dimensions in the plot (profile distribution). 
Usually, as in this case, the two match fairly well. The excep- 
tions arise when there is a high-likelihood region that only 
occupies a small volume of parameter space. 

The two lower panels of Fig[2]show that the three pa- 
rameters controlling the SN feedback are all positively cor- 
related with each other, in a way that allows us to learn 
how the A"-band luminosity function constrains the feed- 
back model. To do this we rewrite equation [4] as 



Am e jectod _ _ V$ N ( \ 6di s k Vyi t 
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= Chalo 
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vir 



7HT > (12) 



£halo y 

from which we see that the amount of ejected gas per unit 
mass of star formation drops to zero for halos with virial 
speed greater than 



Kdr.O = 



( £halo 



Vet 



(13) 



V £disk , 

In our analysis this cutoff is represented by the line of max- 
imum likelihood in the lower left panel of Figure [2] and 
corresponds to ivir ~ 140 km s" 1 , which translates into 
M s tar ~ 10 10 5 M Q and M K » -23. This cutoff virial ve- 
locity is lower than in DLB07, which means that our SN 
feedback stops being effective at fainter magnitudes, allow- 
ing more stars to form in L* galaxies. Since we need to 
assume a stronger SN feedback to decrease the faint end of 
the luminosity function, this is the only way to assure that 
enough stars will still form in brighter galaxies. 
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produced by this process, 



LbH OC /bh &AGN JTlcold /hot • 



(14) 
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Figure 2. Correlations between the 6 parameters analysed in our 
study only constrained by the \ 2 test on the i\-band luminosity 
function. For the values of the parameters plotted in log space, 
the solid contours represent the 68% and 95% preferred regions 
from the MCMC (the posterior distribution) and the colours the 
maximum likelihood value sampled in each bin (the profile dis- 
tribution). The colour scale is normalized by the maximum like- 
lihood value of 0.87. White regions in the plot represent regions 
either with very low likelihood, less than 0.1 per cent of the peak, 
or regions that have not been visited by the MCMC chain. 



For a given value of Vvir,o the amount of ejected gas is 
proportional to thaio, with the maximum likelihood solutions 
showing a linear relation 7 c j oc ehaio- This corresponds to a 
roughly constant amount of gas being held in the external 
reservoir: in a steady state the external gas content is pro- 
portional to the ratio of the influx and outflux rates. Our 
regions of high likelihood represent a considerably higher 
amount of gas being in the external reservoir than in DLB07 
with a corresponding reduction in star formation in faint 
galaxies. 

The value of edisk, controlling the reheating of cold to 
hot gas, has a minor impact except as a way of controlling 
the critical magnitude limit above which feedback is inef- 
fective. Presumably cooling times are so short in galaxies 
below the magnitude limit that any gas that is reheated will 
quickly cool down again. One might have expected that it 
could be used to control the stellar mass fraction in large 
galaxies where the cooling time is relatively long, however 
this does not seem to be the case. 

The AGN feedback parameters, shown in the top-left 
panel of Figure [2j have a broader acceptable region but also 
show a high-likelihood spine that runs diagonally down from 
top-left to bottom right. This can be explained by combining 
equations [6] [7] and [8] to obtain the mechanical heating rate 



Thus, for given cold and hot gas fractions in the galaxies, the 
line represents a single heating rate. This degeneracy is bro- 
ken if the BH masses are used as a constraint (section 14.41) . 
since their growth is mainly dominated by the quasar mode. 

The maximum likelihood channels described above all 
have star-formation efficiencies of esF ~ 0.04, similar to that 
of DLB07. This reinforces the conclusion that the SN and 
AGN feedback parameters act to maintain a constant mass 
of cold gas available for star-formation. 

Apart from the main band discussed above, the two up- 
per panels in Figure [2] show alternative regions of acceptable 
likelihood. These have lower star-formation efficiencies, re- 
quiring a greater mass of cold gas. This, in turn, leads to 
a smaller product of /bh and A?agn- We do not dwell on 
these solutions here as they seem to be ruled out by other 
observations: in particular they result in excessive BH-bulge 
mass ratios. 

4.3 The Colour-Stellar Mass Relation 

The star formation rate is another essential quantity in char- 
acterising the galaxy population. Although this property can 
be directly extracted from the models, it is not easily com- 
parable with observations. While in the models, the mass 
transformed into stars at each time step is computed, only 
indirect observational estimators are available. 

Galaxy colours are one such indirect measure of the re- 
cent star-formation history of galaxies, with a clear bimodal- 
ity between an old, passively evolving re d population, and 
a young star forming blue sequence (e.g. Kauffmann et all 
2003!; iBrinchmann et al I l2004l : iBaldrv et ai]|2004h . 

The top panel of Fig. [3] shows the colour-magnitude 
relation for the DLB07 model. Although it has some prob- 
lems in correctly predicting the slope and fraction of each 
population in some m ass ranges (see IBaldrv et af] 120061 : 
I Weinmann et all 120061 ) . it clearly reproduces the bimodal- 
ity. 

To test the correctness of model colours, we divide the 
galaxi es into the two populat ions using the selection crite- 
ria in [Weinmann et al.l (120061 ). (g - r) = 0.7 - 0.032 (M r - 
5 log ft+16.5), converted into a cut on the colour-stellar mass 
relation at redshift zero, (B — V) = 0.065 log(M*/i 2 / M ) + 
0.09, and shown as the solid line in the upper panel of Fig. [3] 
The conv ersion from the g — r t o the B — V colour was done 
following iFukugita et al] (jl996ft , a-r = 1.05(B - V) -0.23. 
The fraction of red galaxies for dif ferent mass b i ns is t hen 
compared with observations from IBaldrv et af] (|2004T l a.s 
shown in the lower panel. The observatio nal masses based 
on the 'diet' Salpeter IMF i|Bell et al.ll2003l ) were reduced by 
0.15 de x to agree with the IMF assumed in our semi-analytic 
model (|Chabrienl2003h . The fact that a different blue band 
was used in the observational colour cut (it — r instead of 
our converted g — r) could potentially lead to discrepancies 
in the number of objects identified in each population if the 
two gaussian distributions defining each population signif- 
icantly overlap. However, as Fig. [3] shows, there is a clear 
division between red and blue galaxies in the B — V CM di- 
agram for model galaxies. In this way any small differences 
in the number of red galaxies caused by the different colour 
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Figure 3. Comparison between model and observational colours. 
The upper panel shows the DLB07 B-V colour-stellar mass rela- 
tion, with the solid line representing the empirical division be- 
tween populations in the model. On the bottom panel the frac- 
tion of red galaxies as a function of stellar mass from DLB07 (red 
open circles) and our b est fit model (red fille d circles) is compared 
with observations from Baldrv et al ] l|2004) (filled blue squares). 



cut used, should be well within the 0.05 error assumed for 
the model frac tion. 

Following ICroton et alj (|2006l ) we take the resolution 
limit for model colours to be at a stellar mass of approx- 
imately 10 9 ' 5 /i _2 Mq, above which the model reproduces 
the red fraction reasonably well, with some minor excess for 
galaxies above L*. 

The agreement between model and observational 
colours is calculated using a maximum likelihood method 
with a constant value for the errors in the model (cr mo de! = 
0.05) given by the variation in the fraction of red galaxies in 
a sample of 20 sub-volumes of the Millennium Simulation, 
similar to the one used in our analysis. 




Figure 4. As for Figure[2] but constrained only by the maximum 
likelihood test on the fraction of red galaxies. The colour scale is 
normalized by the maximum likelihood value of 0.89. 



We assume that both model and observational values 
are Gaussian distributed around the true fraction F, with a 
likelihood 



C 



(Colour) — exp 



(/model - F) (/ obs - Ff 



2(T 2 

LO model 



that has a maximum value 



-C(Coiour) — exp 



/ob 



2K 



(15) 



(16) 



The DLB07 model, tuned to reproduce a different set 
of observational colours (the red and blue luminosity func- 
tions from 2DFGRS), correctly predicts the fraction of red 
galaxies at low and high masses but makes the transition 
from moderate to high red fractions at a stellar mass which 
is too low. 

In Fig [4] we plot the allowed regions in likelihood and 
posterior space. Perhaps surprisingly, the colour constraint 
picks out a similar relationship between ehaio and edisk as 
does the J\"-band. This is because it also requires a cessation 
of SN heating in galaxies with virial speeds above 140 km s" 1 
which would otherwise have an excessive red fraction. 

The constraint again requires a constant mechanical 
heating from AGN feedback (as shown by the line of maxi- 
mum likelihood in the upper-left panel of the figure) which 
is responsible for the elimination of blue galaxies at high 
masses. The line of highest likelihood lies slightly below that 
seen for the /^-band constraint. Along with this, the upper- 
right panel of the figure shows a preference for a slightly 
lower star-formation efficiency. However, in each case there 
is an acceptable region where the allowed parameter spaces 
overlap. This changes when we move to our third constraint. 
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Figure 5. The BH-bulge mass relation for DLB07 (contours 
and dots) is comp ared with local observational data from 
lHaring fc R ix (2004) (blue crosses). The best fit to the observa- 
tional data points is given by the blue line running from bottom- 
left to top-right, while the two black lines perpendicular to this 
relation divide galaxies into the two mass bins used to compute 
the likelihood. 



4.4 The Black Hole-Bulge Mass Relation 

We have seen in the previous section that the power of the 
radio-mode AGN feedback depends upon the product of the 
quasar and radio-mode growth factors. However, the mass 
growth of the BHs is dominated by the quasar mode alone. 
We can therefore use the BH-bulge mass relation to break 
this degeneracy. 

If we require semi-analytic galaxies to be constrained 
solely by this relation, then the other parameters in the 
model will be free to shift into implausible values allowing 
any point in parameter space to have a reasonable likeli- 
hood. For this reason, we require model galaxies to follow 
both the BH-bulge mass relation and the if-band luminos- 
ity function of bright galaxies (i.e. the host galaxies of these 
BHs). 

The BH and bulge masses for the original model are 
plotted in Fig. |5]and compared with local observations from 
lHaring fc Rb3 ?2004). The model galaxies fall on top of the 
observational best fit (given by the blue line) , with the scat- 
ter in the relation also reproduced. 

In order to test the semi-analytic results against ob- 
servations, we divide the data into two bins (perpendicu- 
lar to the observational best fit), represented by the solid 
black lines on Fig.0 15.2 < ttibh + 0.90 m Bu igc 17.75 and 
17.75 < raBH + 0.90mBui g e, and for each of the bins we com- 
pute the binomial probability for the observed distribution 
of mass ratios above and below the best-fit line, given the 
fractional distribution from the model galaxies: 



-(BH-Bulgc) 



2I p (k,n-k + l), 
2(l-Jp(*,n-* + l)), 



4<0.5 (17) 
J p >0.5 1 ' 
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Figure 6. As for Figure[2] but constrained only by the binomial 
test on the BH-bulge mass relation and by the i\-band luminos- 
ity function of galaxies above M k = —23. The colour scale is 
normalized by the maximum likelihood value of 0.86. 



where k is the number of observed galaxies above the best 
fit in each bin, n the total number of observed galaxies in 



the same bin, and p is the equivalent fraction, k/n, for the 
model galaxies in th e bin. I v (a, b ) is t he incomplete beta 
function as defined in iPress et all |2007l ). The two formulae 
are required since we need to exclude both extremes of the 
distribution, corresponding to an excess of points both above 
and bellow the best-fit line. 

In Fig. [S] we plot the posterior and profile likelihood 
distributions of the parameters constrained only by this bi- 
nomial test and the _K"-band luminosity function of galaxies 
brighter than Mk — —23. 

It is immediately apparent that the region of high like- 
lihood is much smaller for this test than for the first two 
constraints. Moreover, this region corresponds to high val- 
ues for the SN feedback parameters, combined with a low 
AGN feedback efficiency. As such, it is incompatible with 
the acceptable regions in those previous tests. 

However, there is a lower-likelihood, but still acceptable 
(likelihood > 0.1) region that extends towards lower SN and 
higher AGN parameters. As the MCMC contours show, this 
occupies a much larger volume of parameter space than the 
high-likelihood peak (and so, in a Bayesian sense, the true 
solution is more likely to be found in the former than the 
latter) . 

Looking at the upper-left panel in the figure, we can 
see that the acceptable region runs from bottom-left to top- 
right, perpendicular to the lines seen in the previous two 
tests. The BH-bulge mass ratio thus breaks the degeneracy 
in the AGN parameters. 



MCMC in Semi- Analytic Models of Galaxy Formation 9 




Figure 7. As for Figured] but constrained by all three observa- The solid lines represent the maximum likelihood in each bin 

tional properties: the K-band luminosity function, the fraction of marginalised over the other dimensions in parameter space, 

red galaxies, and the BH-bulge mass relation. The colour scale is 
normalized by the maximum likelihood value of 0.037. 



5 COMBINED OBSERVATIONAL 
CONSTRAINTS 

The likelihood of the model for a given point in parameter 
space is computed by taking the product of the three inde- 
pendent observational constraints described in the previous 
sections. 

1i(xi) = £(K-Band) X ^(Colour) x £(BH-Bulgc) (18) 

This quantity is calculated at each MCMC step and 
used to derive the acceptance probability (equation I10[). We 
run our sampling over approximately 30 000 steps, exclud- 
ing an initial burn-in to assure the independence of the final 
results from the starting point of the chain. The output is 
then analysed using getdist, part o f the COSMOMC soft- 
ware package l|Lewis fc Bridle! |2003 >. As for the individual 
constraints, the code is adapted to produce Id and 2d max- 
imum likelihood (profile) and MCMC marginalised (poste- 
rior) distributions, best values and confidence limits for the 
parameters, convergence statistics for the chain and corre- 
lation statistics for each parameter. 

Although, in our sampling we do not impose rigid limits 
on the parameter range, Figure [7] shows that the preferred 
regions are well-constrained. Figure [S] shows the profile dis- 
tribution for each parameter, marginalised over all others. 
Unfortunately, the maximum likelihood is just 0.037, thus 
the best-fit solution is incompatible with the three combined 
observations at the 2-cr level. 

The principal cause of the low likelihood is an incom- 
patibility between the BH-bulge mass constraint and the 
K-haxiA and B ~ V colours. As discussed in section [6] this 
might be caused by observational uncertainties in the black 



hole and bulge masses or might reflect a deficiency in the 
black hole growth model, which is still very simplistic. Nev- 
ertheless, observational uncertainties, principally that asso- 
ciated with stellar population synthesis modelling, make it 
premature to conclude that the DLB07 formalism is ruled 
out. 



5.1 Best Fit Parameters and Confidence Limits 

The best fit and confidence limits for the 6 free parameters, 
together with the published values from DLB07 are shown 
in table [2] All the parameter values in the original model 
fall within our 2a regions except for the SN disk reheating 
efficiency, which we require to be larger than before. 

Both the star formation efficiency and the AGN quasar 
mode parameters from DLB07 closely match our best fit 
values, while our AGN radio mode efficiency is slightly lower 
than before. 

For the SN feedback parameters, the original halo ejec- 
tion efficiency, thaio is below our best fit, whereas the origi- 
nal gas reincorporation efficiency, 7 j , is slightly higher. This 
combination acts so as to produce a higher fraction of gas 
trapped in the external reservoir with the new parameters, 
and hence a smaller mass of cold gas available for star- 
formation in dwarf galaxies. 

The DLB07 SN reheating efficiency tdisk, is consider- 
ably lower than our best- fit. As discussed earlier, the main 
effect of this is to raise the ratio edisk/ehaio with the new pa- 
rameters and hence, from equation 1131 to lower the critical 
virial speed above which feedback is ineffective. 
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Table 2. Statistics from the MCMC parameter estimation for the 6 parameters selected from the original model. The best fit and 
confidence limits (derived from the colour contours in Fig. [7} are compared with the published values from DLB07. 





DLB07 


Mean 


lower (2<r limit) 


lower (lcr limit) 


upper (lcr limit) 


upper (2a limit) 


a SF (SFE) 


0.03 


0.039 


0.020 


0.020 


0.11 


0.13 


kAGN (AGN Radio) 


7.5 x 10~ 6 


5.0 x 10~ 6 


2.4 X 10~ 6 


2.4 X 10~ 6 


9.7 X 10~ 6 


1.1 x 10~ 5 


Jbh (AGN Quasar) 


0.03 


0.032 


0.014 


0.014 


0.103 


0.115 


(disk (SN Reheating) 


3.5 


10.28 


4.43 


4.52 


24.37 


24.37 


(■halo (SN Ejection) 


0.35 


0.53 


0.26 


0.26 


1.17 


1.17 


7ej (SN Reincorporation) 


0.5 


0.42 


0.08 


0.08 


0.73 


0.79 



5.2 Galaxy Properties in our Best Fit Model 

The MCMC parameter estimation was carried out using 
only one data-file representing 1/512 of the Millennium vol- 
ume. In this section we present results using our best-fit 
parameters in the full volume. 

5.2.1 Galaxy Luminosity Functions 

The requirement that star formation be inefficient in low- 
mass galaxies is a common problem to both sets of semi- 
analytic models built upon the Millennium Simulation. 
There is an apparent excess of dwarf galaxie s that can be 
seen in the if -band luminosity function in both lCroton et all 
l|2006l ) and lBower et al l l|2006l ). 

In the left panel of Fig. [5] the 7f-band luminosity func- 
tions from DLB07 and our best-fit model are plotted against 
the observational data set used to constrain the sampling. 
As discussed in section 14.21 to get a good agreement with 
observations, the model needs to form considerably fewer 
stars in low mass galaxies. This is achieved by reducing the 
amount of cold gas available for star-formation by increasing 
the SN heating efficiency and decreasing the amount of gas 
reincorporated at each time step. To assure that enough 
galaxies are produced, the virial velocity cutoff above which 
SN feedback is ineffective is lowered, by raising tdisk relative 
to e h aio (equation [13]). 

In the right panel of Fig(9] the best-fit model seems to 
show poorer agreement with the bj-band than the original 
DLB07. The new fit does reproduce the number-density of 
dwarf galaxies accurately, but shows a large excess around 
L^. This is partly a reflection of the excess seen in the K- 
band in the same region, but has a large magnitude. Given 
the good match to the colour fraction, this seems surpris- 
ing and points to inconsistencies and/or uncertainties in the 
conversion of mass to luminosity via stellar population syn- 
thesis (see section [5~3]) . 

5.2.2 Galaxy Colours 

In Fig. [10] we show the predicted galaxy colours in our best 
fit model, the B-V colour- stellar mass relation on the top 
panel and the fraction of red over the total number of galax- 
ies as a function of stellar mass in the bottom panel. Our 
best fit correctly reproduces the fraction of red galaxies by 
slightly increasing the number of blue galaxies around 
compared to DLB07. 



The colour-stellar mass relation also shows improve- 
ments, keeping the bimodality between the red and the blue 
galaxies, but increasing the slope of each population as sug- 
gested by observations. Nevertheless, near the lower-mass 
limit we impose in our study a population of red dwarfs 
starts to emerge, representing the highest number density 
peak in the red population. This is in disagreement with ob- 
servations, where the majority of the red galaxies are mas- 
sive, and the dwarfs are predominately blue. We address this 
problem, and explore possible solutions in section [S] 



5.2.3 The Black Hole-Bulge Mass Relation 

The best-fit BH-bulge mass relation is almost unchanged 
from that shown in Figure[5]and so we do not repeat it here. 
There is enough freedom in the model to allow the AGN 
parameters to adjust themselves to recover the correct BH 
masses, despite the differences in the SN parameters between 
DLB07 and out best fit. 



5.3 The Galaxy Stellar Mass Function 

As discussed in section |4T^1 the stellar mass function is one 
of the most fundamental properties of a galaxy population, 
but it is difficult to derive accurately from observations. In 
Fig. [TJJwe show how our best-fit model and DLB07 masses 
comp are with observa tion ally-derived stellar m ass functions 
from lBell et~ai1 (|2003h and lBaldrv et al.1 (120081 ). The latter is 
one of the most robust mass derivations, using the New York 
University Value-Added Galaxy Catalogue that combines 4 
different methods for determining galaxy masses from SDSS 
data. The error bars in the figure span the maximum and 
minimum mass estimates from that analysis. 

Aft er all the data s ets are converted into the same IMF 
(that of IChabrieri 12003! ) . the compariso n between our b est- 
fit model and the masses derived from iBell et all (|2003t ) m 
Fig. Illl shows the same behaviour as the Jf-band luminosity 
function. With our more effective supernova feedback the ex- 
cess of dwarf galaxies largely disappears, but there is a slight 
excess of L* galaxies. Small differences might arise when 
comparing both the fc-band and stellar mass functions from 
the model and a specific observational data set, even though 
similar stellar population synthesis models were used. This is 
because, to convert stellar mass into luminosity, one should 
have knowledge about the age and metallicity of the galaxy 
stellar population. These quantities are directly available in 
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Figure 9. Comparison of the predicted _R"-band (left panel) and bj-band (right panel) luminosity functions at z = from DLB07 
(dashed red line) and our best fit model (solid red line). On the left panel, the data p oints represent the observations us ed to constrain 
the luminosities of galaxies in our MCMC parameters estimation llCole et al.ll200ll ; iBell et alj|2003: 1 Jones et al.ll2006l) . On the right 
panel, the 6j-band luminosity function i s co mpared with obser vations from 2DFGRS (green filled squares) and 6DFGS (blue open 
circles), respectively iNorberg et al.l J2002f) and I Jones et al.l d2006T) ■ 



the model, but in observations they are difficult to derive 
and are subject to large uncertainties. 

Also shown in Fig. 1 11 | is a comparison o f the predicted 
masses with data from iBaldrv et all |2008). A systematic 
difference between the model and the data is evident, with 
the former predicting a much larger number of galaxies on 
and above L*. While the hor izontal error bars plotted in 
Fig.fTTIfor lBaldrv et all (|2008h reflect only the bin size, the 
authors refer to differences as large as 0.15 dex in mass es- 
timates from the d ifferent methods. Another recent work 
|Conrov et alj [2008) points to even larger errors of up to 
0.3 dex that may result from imprecise modelling of key 
phases of stellar evolution. If these uncertainties lead to the 
observationally-derived masses of large galaxies being un- 
derestimated b y about 0.15 dex, th en our best fit and ob- 
servations from IBaldrv et al.l |2008h would be in extremely 
good agreement throughout the whole mass range. 

The differences between data sets highlight the need for 
caution when galaxy formation models are compared with 
observations. In principle one should expect the properties 
and allowed parameter ranges to change if either the stellar 
population synthesis or the dust model need to be read- 
justed. In this paper we have chosen to fix these so as to 
focus our study on the parameters controlling the most ba- 
sic properties of the semi-analytic model: star formation and 
feedback. 



6 DISCUSSION AND CONCLUSIONS 

Since they were in troduced as a new technique to understand 
galaxy formation (IWhite fc Reedll97i IWhite fc Frenkll99ll ; 
lLacev fc Silklll99ll ; ICold Il99lh . semi-analytic models have 
always lacked a proper statistical analysis of the allowed 



range of their free parameters and a consistent way to test 
the goodness of the fits produced. 

To overcome this weakness we have implemented a 
Monte Carlo Markov Chain parame ter estimation technique 
into the lDe Lucia fc Blaizotl (|2007T) semi-analytic model, to 
obtain the best values and confidence limits for the 6 free 
parameters in the model responsible for shaping the stel- 
lar mass function and for the colours of galaxies. Com- 
paring the model with three different observational con- 
straints s eparately: t h e combined if -band lum i nosity func- 
tion f rom ICole et all d200lh: IBell et all il200St) ; Ijones et all 
l|2006l ), galaxy colours fro mlBaldrv et all (12004 ) and the BH- 
bulge mass relation from lHaring fc RbJ (|2004| ). we are able 
to identify which particular parameters (and hence which 
galaxy formation processes) are responsible for each individ- 
ual property and which show correlations and degeneracies. 

Combining the three observational tests, we are able 
to fully constrain the model parameters, obtaining a best 
fit and confidence limits within the very limited region of 
acceptable likelihood. Our best model is given by: asF ~ 
0.039lH?9. k AG N = (5.0±1;°) x HT 6 , f BH = 0.032±°;°;i, 



C d / s k 



10.28 



+ 14.09 

5.85 ' t halo 



^+ - 64 
u.oo_ 27 



and 7 e j = 0.42 



+0.37\ 



As shown in table[2] all the parameters in the original model, 
except the supernova reheating efficiency fall within our 2- 
a confidence limits. Our best fit maintains the values for 
the star formation efficiency and for the AGN quasar mode 
parameters, while increasing the SN gas reheating and ejec- 
tion and decreasing the AGN radio mode and gas reincor- 
poration efficiencies. For our preferred set of parameters 
the model has a likelihood of n(xi) — L( Mass ) x L(c iour) x 
i(BH-Buige) = 0.037. This value means that the best-fit so- 
lution is incompatible with the three combined observations 
at the 2- cr level. 

In this paper, we have used the Millennium Simulation 
which adopts a ACDM cosmology. It is possible that the 
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Figure 10. The top panel shows the B-V colour-stellar mass 
relation for the galaxies in our best fit model. The solid line 
represents the division between the red and blue populations in 
IWeinmann et alj [1200611 . The predicted fraction of red galaxies as 
a function of stellar mass is showed in the bottom panel. The orig- 
inal DLB07 model (open red circles) is compared wi th our best fit 
model (filled red circles) and observational data from lBaldrv et all 
||2004T ) (filled blue squares). 



Universe may be better described by an alternative cosmol- 
ogy with fewer low-mass halos. However, our purpose here 
is to try to find an astrophysical solution that is compatible 
with ACDM. 

As discussed in previous chapters, observational uncer- 
tainties, principally that associated with stellar population 
synthesis and dust modelling, make it premature to conclude 
that the DLB07 formalism is ruled out: shifting the observed 
black-hole/bulge mass ratio by 0.15 dex raises the best-fit 
likelihood to 0.07, which is marginally acceptable. Neverthe- 
less, the apparent incompatibility between the black hole- 
bulge mass relation and the other constraints indicates that 



Figure 11. Comparison for the predicted stellar mass function at 
z=0 from DLB07 (dashed red li ne) and our best-fit model (solid 
red line) wi th observations from Bald ry et al. | ll2008h (blue open 
circles) and lBell etldl l|2003l )f green filled squares). 



the black hole growth treatment in the model might be 
too simplistic, in particular by assuming that there is no 
feedback from the quasar mode, when it seems to be re- 
quired to reproduce the X-ray luminosity function of halos 
jBower et al.ll2008l ; [Short fc ThomasfeOOSl ). As discussed be- 
low, some additional recipes might also need to be included 
for the model to better reproduce observational luminosities 
and colours, which could in principle increase the likelihood 
of the best fit model. 

We produce a isT-band luminosity function for our 
best parameter values that improves the agreement with 
observations at the low-luminosity end. This is achieved by 
taking a higher heating efficiency from SN and a lower rein- 
corporation rate of gas ejected from the halo. This reduces 
the amount of cold gas available to form stars, avoiding 
the excess of faint galaxies in the original recipe. More 
effective SN feedback has bee n proposed in the past. For 
example, iBertone et al.1 (|2007l ) studied a wind model that 
improved the number density of dwarfs for both the mass 
and luminosity function while improving the distribution of 
metals. However, the high value of ejection that is required 
by our model seems to be unrealistic when compared with 
observations ljMartinlll999l ). This indicates that additional 
processes such as disruption of sat ellites through tida l 
effects might need t o be included dBullock et al.l l200ll; 
Taylor fc Babull l200ll; iBenson et~all |2002|; I Monaco et all 



20061; 



2006 



Murante et al. 2007 



.Wei nmann e t al. 
Henriques et al.ll2008l ; ISomerville et al.ll2008T l 

For both the luminosity-function and colour con- 
straints, the SN reheating and ejection parameters are 
strongly correlated, which we interpret as an upper virial 
velocity limit of 140 km s -1 for galaxies that can eject mass 
via SN heating. Since we need to assume a stronger SN feed- 
back to decrease the faint end of the if-band luminosity 
function, this relatively low value assures that our SN feed- 
back stops being effective for galaxies with masses above 
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A/* ~ 1O 1O ' 5 M0. This is the only way to ensure that enough 
stars will form in brighter galaxies, and also produces more 
blue galaxies around L t than the original DLB07 parame- 
ters. 

Significant correlations exist between the three param- 
eters governing SN feedback , tdisk, ^haio & 7ej, even in the 
combined analysis. This suggests that the model could be 
rewritten with one or two fewer free parameters. However 
this degeneracy could in principle be broken if the metallic- 
ity of gas and stars were to be considered. 

The model with the original DLB07 parameters cor- 
rectly predicts the bimodality in the colour-stellar mass re- 
lation, however, it has difficulties in matching the exact 
number of the blue and red sequence galaxies. Our best- 
fit model correctly predicts the relative fraction of galaxies 
in each colour population, however the early cutoff on SN 
feedback leads to an excess of galaxies with masses between 
10 10 5 and 10 110 M Q . Furthermore, as the original model, it 
shows a large population of small red galaxies in the B — V 
colour-stellar mass plot in contradiction with observations. 

The problems with low-mass galaxy colours in semi- 
analytic models have been identified in the pas t , particularly 
the e xcess of red dwarfs l|Croton et al ] |2006l ; iBaldrv et al.l 
2006). Possible solutions might include the delayed stripping 
of gas from satellites after their dark matter halo is disrupted 
(allow ing them to coo l gas, form stars and stay blue for 
longer lFont et al1l2008l ). Or again, tidal disruption of dwarfs, 
which would affect mostly red, satellite galaxies. This would 
move them to even lower masses, below 1O 9 '°M0 (where an 
upturn in the stell ar mass function is se e n) and produce 
intra - cluster light dWeinmann et al. 2006; [Henriaues et al.l 
l20Qg| ; ISomerville et al.ll2008l ). 

The purpose of this paper is to show that MCMC pa- 
rameter estimation techniques, adapted from those used in 
cosmology, can be used to map out likelihood contours in 
the parameter space of semi-analytic models of galaxy for- 
m ation. For this particular a nalysis we chose the formalism 
oflDe" Lucia fc Blaizotl (|2007l ). but the same method could 
equally be applied to other models. 

In the future, we would like to extend the method to 
undertake model selection, providing an objective measure 
of the relative value of models with different numbers of free 
parameters. 
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